function [equalorder, ...
          Ugridreduced, sgreduced,...
          id_Ugridreduced,...
          s_id_gridreduced, sU_id_gridreduced,...
          n_U_sets, Tcoord, indices_pairs,...
          numeqconstraints,coordrowseq, fillAeq]=useful_anyparam4b(u0grid, u1grid, u2grid, sg)     
%% U sets
n_U_sets=2; %number of U-sets (number of taste shocks - 1)
Ugrid=cell(1,n_U_sets); %cell 1xn_U_sets

Ugrid{1}=[u0grid-u1grid   u0grid-u2grid -u0grid+u1grid ...
         -u0grid+u2grid  ...
          zeros(sg,1) ...
          Inf*ones(sg,1) -Inf*ones(sg,1)];

Ugrid{2}=[u2grid-u1grid -u2grid+u1grid u2grid-u0grid ...
         -u2grid+u0grid...
          zeros(sg,1) ...
          Inf*ones(sg,1) -Inf*ones(sg,1)];


%% Reduce number of grid points to check by applying our partitioning arguments
%equalorder (sgx1) assigns the same index to parameters vectors such that Ugrid{1}, Ugrid{2}, ..., Ugrid{n_U_sets} have the same order

d=cell(n_U_sets,1);
index=cell(n_U_sets,1);

for j=1:n_U_sets 
    [sortedUjgrid,index{j}] = sort(Ugrid{j},2); 
    %sortedUjgrid gives Ugrid{j} with each row sorted from smallest to largest; it is sgxn_columns_Ujgrid
    %for each row i of sortedUjgrid and for k=1,2,...,n_columns_Ujgrid, index{j}(i,k) gives the column position in Ugrid{j}(i,:) of sortedUjgrid(i,k); it is sgxn_columns_Ujgrid
    d{j} = diff(sortedUjgrid,[],2) == 0; %sg x (n_columns_Ujgrid-1)
    %Suppose that Ugrid{j} has 3 columns. Then, for each row i of sortedUjgrid, d{j} will be: 
    %[1 1] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[1 0] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
    %[0 1] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[0 0] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
end

[~,~,equalorder] = unique([[index{:}] [d{:}]],'rows', 'stable'); %sgx1 

%Reduced grid of parameter vectors and U sets
[~,index_final,~]=unique(equalorder, 'stable'); %one row from equalorder per ordering group together with its index
Ugridreduced=cell(1,n_U_sets);
for j=1:n_U_sets 
    Ugridreduced{j}=Ugrid{j}(index_final,:); %sgreduced x n_columns_Ujgrid
end
sgreduced=size(Ugridreduced{j},1); %size of the reduced grid of parameter values


%% Construct indices to keep track of equal elements inside Ugridreduced{j} for j=1,...,n_U_sets
id_Ugridreduced=cell(n_U_sets,1);
for j=1:n_U_sets 
    id_Ugridreduced{j} = NaN(size(Ugridreduced{j})); % %sgreduced x n_columns_Ujgrid
    for k = 1:sgreduced
        [~, ~, id_Ugridreduced{j}(k,:)] = unique(Ugridreduced{j}(k,:), 'stable'); 
    end
end

s_id_gridreduced=zeros(sgreduced, n_U_sets); %useful to set number of unknowns of the linear programming
for j=1:n_U_sets
    s_id_gridreduced(:,j)=max(id_Ugridreduced{j},[],2); %number of distinct elements of Ugridreduced{j}
end
sU_id_gridreduced=prod(s_id_gridreduced,2); %sgreducedx1 reporting the cardinality of the Cartesian product of the U sets
                                            %sU_id_gridreduced(g) is the number of unknowns of the LP for each g-th parameter value
s_id_gridreduced_2=zeros(1, n_U_sets); 
for j=1:n_U_sets
    s_id_gridreduced_2(:,j)=size(id_Ugridreduced{j},2); %number of elements of Ugridreduced{j} (count as separate even equal elements)
end
sU_id_gridreduced_2=prod(s_id_gridreduced_2); %cardinality of the Cartesian product of the U sets   
%% Invariant components of the linear programming problem 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vectors = arrayfun(@(x) {1:x}, s_id_gridreduced_2); 
n = numel(vectors);
Tcoord_temp = cell(1,n);
[Tcoord_temp{:}] = ndgrid(vectors{:}); 
Tcoord_temp = cat(n+1, Tcoord_temp{:});
Tcoord = reshape(Tcoord_temp,[],n); %matrix of dimension sU_id_gridreduced_2xn_U_sets reporting the column coordinates of all the possible n_U_sets-tuples from the U-sets
%For example if n_U_sets=3 and s_id_gridreduced_2 is [2 3 1], then 
%[ca, cb, cc] = ndgrid(1:2, 1:3, 1:1);
%Tcoord = [ca(:), cb(:), cc(:)];   
%Tcoord= [1     1     1
%         2     1     1
%         1     2     1
%         2     2     1
%         1     3     1
%         2     3     1];
indices_pairs=pairIndices(sU_id_gridreduced_2); %row indices of all possible unordered pairs of rows from the matrix obtained by taking the Cartesian product of the U sets
                                         %[sU_id_gridreduced_2*(sU_id_gridreduced_2-1)/2]x[2]   
%Continuing the example above
%indices_pairs=[1 1; 1 2; ...; 1 27; 2 3; ...; 2 27; ... 26 27]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numeqconstraints=3+14+6;

coordrowseq=[1 1 1 ... %observational equivalence
             2 2 2 ...
             3 3 3  ...
             4:1:3+13 ... %zeros
             3+13+1 ... %ones 
             3+13+1+1:1:3+13+1+2 ... %symmetry
             kron(((3+13+3)+1:(3+13+3)+4), ones(1,2))]; 
              
         
fillAeq=[1 -1 -1 ... %observational equivalence
         1 -1 -1 ...
         1 -1 -1 ...
         ones(1,13) ... %zeros
         1 ...    %ones
         1 1 ... %symmetry 
         repmat([1 1],1,4)];
     
end